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We calculate the normalized second-order correlation function for a system of two tunnel-coupled 
photonic resonators, each one exhibiting a single-photon nonlinearity of the Kerr type. We employ 
a full quantum formulation: the master equation for the model, which takes into account both 
a coherent continuous drive and radiative as well as non-radiative dissipation channels, is solved 
analytically in steady state through a perturbative approach, and the results are compared to exact 
numerical simulations. The degree of second-order coherence displays values between and 1, 
and divides the diagram identified by the two energy scales of the system - the tunneling and the 
nonlinear Kerr interaction - into two distinct regions separated by a crossover. When the tunneling 
term dominates over the nonlinear one, the system state is delocalized over both cavities and the 
emitted light is coherent. In the opposite limit, photon blockade sets in and the system shows an 
insulator-like state with photons locked on each cavity, identified by antibunching of emitted light. 

PACS numbers: 42.50.Ar, 42.50.Pq, 71.36.-|-c, 73.43.Nq 



I. INTRODUCTION 

Recent advances in cavity quantum electrodynamics 
(CQED) have led to the demonstration of a number of 
striking phenomena related to the fundamental proper- 
ties of light-matter coupling, e.g. when single or few 
quantum emitters interact with the mode of an elec- 
tromagnetic resonator!^*'^' The experimental realization 
of the Jaynes-Cummings (JC) model;'' which predicts a 
strong light-matter coupling regime when the Rabi fre- 
quency between the oscillators exceeds their respective 
loss rates;^ ^ has been achieved in both atomic^ and solid- 
state^ CQED with single two-level emitters in high-Q 
resonators. In the strong coupling regime, the CQED 
system is intrinsically anharmonic at the level of single 
quanta, which derives from the underlying anharmonic 
nature of the emitter's eigenstates.lSHSl 

Following these early works, the ultimate limit of non- 
linear optics, i.e. the ability to control the nonlinear 
response of a system by the injection of single photons 
through the so called photon blockade effect,!^ has been 
recently reached. Inhibition of the resonant transmission 
of a single photon because of the presence of another 
one within the cavity has been experimentally demon- 
strated with both single atoms in optical cavitie&iS and 
semiconductor quantum dots strongly coupled to pho- 
tonic nanocavities.'i^l In all these experiments, the statis- 
tical properties of the resonant light beam transmitted 
through the nonlinear system gives precise information 
on the nature of the effective photon-photon interaction 
within the cavity. Photon blockade has been shown to 



be strictly characterized by conversion of a classical, co- 
herent field at the input into a nonclassical, antibunched 
photon stream at its output. ^^ Similar effects have been 
also predicted for other types of coherent fields .'^^"^Sl 
More recently, there has been an intense effort towards 
the exploitation of nonlinearities arising from Coulomb 
interaction in confined electron and photon systems.'^ In 
such a case, mixed light-matter states (polaritons) arise 
from the strong coupling regime of quantum well excitons 
and microcavity photons, which are three-dimensionally 
confined thanks to the progress in lithographic tech- 
niques. The polariton quantum blockade has been pre- 
dicted for such highly nonlinear light-matter stateSjE^and 
first evidences of nonlinear behavior have been reported 
experimentally.'22l These systems are likely to provide a 
further playground for single-photon nonlinear optics in 
the near future. 

Motivated by the great level of control achieved in 
CQED experiments with single cavities, recent theoret- 
ical work has explored multi-cavity nonlinear systems. 
Initial work was mainly aimed at stu dying t he superfluid- 
insulator quantum phase transitiorP^HSSl of the effec- 
tive Bose-Hubbard hamiltoniaUjlSZHlII or the JC-Hubbard 
model22tt29] for arrays of CQED systems under quasi- 
equilibrium conditions. Subsequent work dealing with 
coupled non-linear cavity syste ms ha s addressed the dy- 
namics in a two-site JC model,H2IS] soliton physics,H2l a 
proposal for observing fractional quantum Hall states,l331 
the possible realization of a Tonk s-Gir ardeau gas in dif- 
ferent one-dime nsiona l geometries ,'2^'^ the study of effec- 
tive spin modelaSSEZl and of entanglement generation,SSl 



the use of coupled cavity systems as efficient single- 
photon sources even in the presence of weak pho- 
ton nonlinearitiespS and the signatures of superfluid- 
insulator quantum phase transition for an infinite CQED 
array under pulsed coherent driving. ^"^ 

In a recent work,'^^ a proposal has been made to ob- 
serve signatures of strong photon correlations in a sys- 
tem of three coupled nonlinear cavities, with the cen- 
tral one displaying single-photon nonlinearity, through 
the measurement of its degree of second order coher- 
ence. Besides being a readily realizable system with 
state-of-the art techn ology with both atomic and solid 
state CQEDpEillllin] ^j^jg gygtem is a possible quantum 
photonic device in which information encoded in classical 
field states can be controlled by the presence or absence 
of single photon quanta (a single-photon transistor) . 

In order to extend and generalize the latter work, here 
we present a systematic theoretical analysis of a model of 
two coupled cavities, both of them assumed to be nonlin- 
ear at the single photon level. The model takes into ac- 
count coherent driving as well as global dissipation chan- 
nels within a master equation treatment (see a sketch of 
the system, Fig. nj. The dynamical equilibrium reached 
by the system is due to the balance between pumped and 
dissipated photons and it is analyzed in steady state. Fol- 
lowing previous work,'^ we mainly concentrate on calcu- 
lating the second-order correlation function for the light 
emitted from each cavity, both analytically and numeri- 
cally. Interplay of coherent tunneling and on-site interac- 
tions is clearly identified in the crossover from Poissonian 
to sub-Poissonian light statistics of the emitted light. 

The paper is organized as follows: we first introduce 
the model and the master equation in Sec. II. In Sec. Ill 
we provide a description of the analytical solution for 
the master equation that can capture the second-order 
correlation function for this model, and compare it to a 
full numerical solution. 



II. 



THEORETICAL FRAMEWORK 



We will investigate here the photon correlations in the 
two-site Kerr-Hubbard model (Fig. [I]) given by {h — 1) 



Hkh = Y. ^^^pIp^ + U^p\p\p^P^ + F^e-'^^'p] + F*e'^^'p, 

i=l,2 

+ j{p\p2+pIpi), (i; 



where Pi {pi ) destroy (create) generic bosonic excitations 
in each of the two cavities at their fundamental frequen- 
cies uJi, Ui is the nonlinear Kerr-type interaction in each 
cavity, J is the inter-cavity tunneling rate, and Fie~^'^^^ 
is the coherent driving amplitude at laser frequency wl. 
There are several ways to realize the model in Eq. ([I]) , 
with either atomic or solid-state cavity QED technology. 
All these realizations principally rely on the formation 
of well-defined quasi-particles, polaritons, of mixed light- 
matter nature, for the availability of the required effective 




Figure 1: (Color online) The system of two coupled nonlin- 
ear cavities. The relevant parameters of the model are indi- 
cated, namely the coherent pumping (F) and dissipation (7p) 
rates, respectively. The cavities are supposed to be nonlinear, 
with an interaction energy U, and they are tunnel-coupled 
by evanescent overlap of their cavity modes with a coupling 
constant J. We assume symmetric cavity parameters in this 
work. 



quasi-particle interactions Ui. One possible way is, e.g., 
to start from two atoms or quantum dots strongly cou- 
pled to their respective cavity modes, with the two cavi- 
ties in optical contact with each other. While the single- 
photon nonlinearity derives from the light-matter cou- 
pling in the framework of a JC model, the tunnel coupling 
is due to photon tunneling and the coherent pump acts di- 
rectly on the photonic degrees of freedom. In the absence 
of losses or spontaneous emission decay, such a system 
would be described by a two-site JC-Hubbar d Ham ilto- 
nian, extensively discussed in recent literature !'^^*^^" ^ In 
the weak pumping limit and the dispersive regime, the 
JC nonlinearity can be reduced to an effective Kerr-type 
nonlinearity between polaritons.'^S 

A more straightforward and conceptually simple way 
of obtaining a Kerr-type nonlinearity that is effective 
at the single photon/polariton level is to consider solid- 
state systems in which the Coulomb interaction is strong 
enough. In particular, we refer here to excitons in quan- 
tum wells coupled to a single photonic mode of a micro- 
resonator where excitons interact via their dipole field.'^S 
In such a case, the Hamiltonian for the two-site system 
in the rotating wave and electric dipole approximations 
is given by H — Hq + Hi witfPS' 

Ha = Yl Kav,»a|a, + uj^,iX\Xi + n,{a\x, + a,x\)] (2) 

i=l,2 

Hi = J2 [V^Xjxjx,X, + E,{t)e-"^^'a\ -f S*(i)e'"^*a,] 

i=l,2 

+ j{a\a2 + a\ai). (3) 

Here, a\ (a^) creates (destroys) a photon in cavity i at 
frequency Wcav,i; while the operators Xi {Xl) describe 
excitonic quasi-particles with energy Wx.i- We assume 



an interaction energy Vi between excitons deriving from 
a contact-type Coulomb interaction, and the exciton- 
photon interaction strength is given by the Rabi fre- 
quency J7,;. Cavity photons are coherently pumped into 
each cavity with amplitudes £^^6"'"""*, and j is the tunnel- 
ing amplitude for photons between neighboring cavities. 
With respect to the latter model, polaritonic excita- 
tions can be defined as linear combination of excitons 
and cavity photons as 



P+.i J \V^ Ui 

where the coefficients are^ 
1 

; Vi = 






(4) 



This transformation then diagonalizes Hq 

Ho=J2Yl ^^AKd^P^.^ (6) 

a=±i=l,2 

with the lower and upper polariton energies respectively 
given by w±,i = (wcav,i+CL;x,i)/2± V^f -I- (Ai/2)2, where 
Ai = Wcav.i — i^x,!- Assuming quasi-resonant pumping of 
the lower polariton level (wl ~ ^-)i Aj < and ne- 
glecting non-resonant contributions, the resulting effec- 
tive Hamiltonian H can be written in the form of Eq. (Ilj) 
with pi = P^^i, iOi = uj^j, Fi = ViEi, J = viV2J and 
Ui = ufVi- Throughout this work we will be assuming 
identical sites for simplicity and drop the site indices on 
the parameters of the system. 

Losses resulting, e.g., from spontaneous emission decay 
of the excitons or cavity photon leakage can be taken 
into account within the quantum Master equation in the 
Born-Markov approximation for the density matrix of the 
quasi-particles in the system, which is expressed in the 
usual Lindblad form 



d 

— p = i[p,iJKH] +'C(p) 



where i?KH is the Kerr-Hubbard Hamiltonian written in 
the rotating frame with respect to the pump frequency 
resonant with the lower polariton mode (a;_ — wl = 0) 

HkH = R(i)^KHR^(t) = Y. [UP\P\P^P^+Fp\+F*p,] 



1=1,2 



St^ 



J{p\P2+P2Pl), 



(8) 



with R(t) = exp{ia;L(pjpi +'P2P2)i}- The Liouvillian can 
be expressed in the usual Lindblad fornPSl 



^=?E[2^' 



PPi 



p\piP - ppIpi] , 



(9) 



where 7p is the polariton dissipation rate in each cavity.!^ 
We will discuss possible realistic implementation of this 
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Figure 2: (Color online) Schematic energy level diagram and 
rates for the coupled cavity system. We use the following 
shorthand notation for the eigenstates: jO,0) — >■ |1), |1,0) — >■ 
12), jO, 1)^13), |2,0)^|4), |1, 1)^15), |0,2)^|6). 



III. RESULTS 



The aim of the present work is to assess the second- 
order correlation function as a quantitative probe of the 
interplay between tunneling and interactions in a two- 
site CQED system described by Eq. (l7|. To this end, we 
calculate the steady state normalized degree of second- 
order coherence,'^ i.e. gL (t), for the light emitted from 
the system, which is defined as 



9ilHr)-9''^ 



{t — ?> GO, r) 



{p^t)p\t + T)p{t + T)p{t)) 

(pHmt))' 



(10) 
(7) In the following, we will only be concerned with the zero- 



model and the required tolerances in Sec. IV 



time delay correlation function in steady state, gss {0) 
{p^"^ p^) / {p^ p)'^ ■ This correlation function of the polari- 
tons can be straightforwardly related to the correlation 
function of cavity photons via Eqs. Q, which is t he q uan- 
tity ultimately measured in a typical experiment .ISIl 

In the weak pumping limit, an analytic solution to the 
steady state master equation for our model can be found 
in the following way. We write Eqs. (l7| and ^ in the 
Fock basis {|ni,n2)}, where rii and n2 indicate polariton 
occupations in cavities 1 and 2, respectively. We consider 
the low energy excitations of the Hamiltonian (Is]) , Ntot — 
ni + n2 < 2. The corresponding energy level diagram and 
rates are schematically shown in Fig. [2j The resulting 
equations of motion for the 36 elements of the density 
matrix can be solved using perturbation theory and a 
recursive procedure in F/jp, as described in App. A. 

The steady-state second-order correlation function 



(2) 

gls (0) for the ith cavity can be calculated as 



(2)^n\ ^ Trj^pIpiKPss} _ E 



5^(0) 



imn' rrn,m' 



Mp\p\piP^W) 



Tr{plp,psJ]2 



E 



st^ 



ra.m' rm,m' 



m\p]p,\m') 



(11) 

where \m) = \n1n2) is a collective notation for the eigen- 
states of the coupled cavity system and p^ ^, is the 
steady-state density matrix calculated from Eqs. Q. We 
will henceforth drop the subscript "ss". In the weak 
pumping limit F/jp <C 1 



-,(2)/ 



,(2) 



<?r (0)= 3^(0) =5(^^(0) - 






(12) 



The calculations are lengthy but straightforward (see 
App. A), and the analytic expressions for the matrix el- 
ements (04^4 and p2,2 are given in Eqs. (A49) and (A29l, 
respectively. Notice that the explicit analytic expression 
for p4_4 is recursively obtained through the analytic ex- 
pression for the elements in Eqs. (A29I, (A33), (A34|, 



(A41), and (A42). The procedure can be generalized, 
e.g., to perturbatively calculate ^^^•'(O) for generic multi- 
site CQED systems, either analytically or numerically 
through the implementation of this recursive algorithm. 
Wc notice that setting J = in the above elements, Eq. 
(12) gives the correct limit of the single cavity photon 



blockade 



.g'^Ho) 



1 + 4([//7p)2 



(13) 



For J 7^ 0, the resulting behavior of ^'^-'(O) is shown 
in Fig. ^ as a function of dimensionless quantities C//7p 
and J/7p. The most striking feature of this plot is the 
sharp boundary that divides regions where g'^'^\0) ~ 
and (7^^''(0) ~ 1. As clearly seen in Fig.pl when the tun- 
neling term dominates over the on-site interaction en- 
ergy the emitted light is Poissonian, ^^^■'(O) ~ 1. This 
reflects the statistics of the coherent driving fields im- 
posed on the system due to the dominance of the linear 
tunneling terms over the non-linear interaction terms in 
the Hamiltonian. Hence the system state is coherently 
delocalized over both cavities with symmetric and anti- 
symmetric combinations of the bare polariton modes. In 
the opposite limit, the emitted photons are antibunched, 
(;(2)('Q~) r^ Q^ which is a clear indication that the number 
of quasi-particles in each cavity can only fiuctuate be- 
tween zero and one. The latter is the photon-blockade 
regime that would be present for J — Q (i.e. individual 
cavities) for U :^ ^p (See Eq. ^^). Note that for J 7^ 
the cross-over takes place for larger values of [/ as J is 
increased, thus showing the interplay of tunneling and 
interactions in the steady state. The boundary of the 
crossover in U, Ui,{J), can be clearly identified following 
the g^^^O) = 0.5 contour (see red dashed fine in the plot) 
and is found to increase approximately linearly with J for 
J,U > -fp. 

In order to check the accuracy of the analytical solution 
presented above, we have numerically solved Eqs. (|7])-(|9]). 
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Figure 3: (Color online) The second-order correlation function 
for the two coupled nonlinear cavities, g^^' (0), as a function of 



U/jp and J/fp (a) as derived analytically from Eq. ( 12 1 and 
(b) from numerical solution of Eq. u\ calculated for F/'fp = 
0.1. The dashed lines in both (a) and (b) show the boundary 
curve Ub{J) defined by the condition 2p4A/p2,2 = 0.5. 



The second order correlation function is numerically cal- 
culated using Eq. (11), by using up to 6 photons in the 



Fock basis to ensure full convergence (a small F/jp is 
assumed to compare with the perturbative analytic solu- 
tion) . The result is shown in Fig. IsFb) with a color scale 
plot, displaying 5^^^(0) as a function of U/^p and J/^p 
to be directly compared to Fig. Isla). To better show the 
agreement and the accuracy of the analytical procedure 
provided in this work, we give in Fig. HI several cuts of 
the color plots of Fig. |3] The quantitative behavior of 
5(^^(0) is very well reproduced by the analytic solution 
over several decades considered in the parameter space. 



IV. PHYSICAL REALIZATION 

The Kerr-Hubbard Hamiltonian Eq. (IT]) can be realized 
with a number of systems, including strongly-coupled 
atom-cavity systems, ^^ circuit QED systems,"^'' or quan- 
tum dots coupled to semiconductor resonators.i^ 
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laritons) with contact-type interaction of Coulomb na- 
ture. In the first case, sizeable values of the nonlinear- 
ity have been al ready shown experimentally in solid-state 
systemspiESE^land tunability might be achieved through 
external control of the qubit-cavity detuning (see, e.g.. 
Supplementary Information Section of Ref. I5ip . In the 
second case, for which photonic crystals represent an 
ideal platform for a prospective implementatiouj-^^ very 
large values of exciton-exciton interactions in confined 
systems have been predicted for single resonators ,'22 and 
experimentally measured recently.ESJ 

Finally, detection of emitted radiation and subse- 
quent measurement of correlation functions would be 
performed by standard Hanbury-Brown-Twiss technique, 
details depending on the geometry of the system un- 
der consideration. For example, in the case of pho- 
tonic crystal slab cavity systems emitted light might 
either be guided through access waveguides and then 
detected in an "in-plane" geometry, or directly im- 
aged in a camera above the slab as in usual micro- 
photoluminescence experiments. Recently, improved far- 
field collection from high-Q photonic crystal cavity modes 
has been reported,!^ which might also be crucial for the 
present proposal. 



CONCLUSIONS 



Figure 4: (Color online) A direct comparison of numerical and 
analytic results for (;"^(0): (a) fixed U/jp and varying J/jp, 
and (b) fixed J/jp and varying U/jp. The full lines repre- 
sent the analytic solutions, while numerical ones are shown 
by symbols. 



Here, we discuss possible realization with 
semiconductor-based optical microcavities. Efficient 
coupling of optical modes in photonic resonators in 
evanescent contact with each other has been already 
demonstrated in a number of different systems, includ- 
ing mic ro-pil lars,!^ micro-disks,™ and photonic crystal 
cavities.lSSESl For the latter, the constant improvement 
of post-fabrication techniques t hat a llow to determin- 
istically tune the cavity modea^I'^ as well as their 
coupling to realize the regime of symmetric sites studied 
here. Note that the tolerance in the parameters for 
observing qualitatively the phenomena studied here 
in the symmetric limit is given by the loss rate "fp. 
A wide range of tunability has been experimentally 
demonstrated in such systems, e.g. for tunnel coupling 
rates ( J ^ — 1 meV),'^and Q-factor engineering allows 
cavity loss rates to be limited to few fieVP^ 

Concerning the nonlinearity, as we have seen in Sec.[ll| 
the effective Kerr-Hubbard model can be realized either 
with strongly coupled two-level systems (qubits) and cav- 
ity modes, or through coupled light-matter states (po- 



The two-site Kerr-Hubbard model represents a 
paradigmatic system to investigate the interplay of tun- 
neling and strong correlations in a CQED system driven 
out-of-equilibrium. We have analyzed this system under 
coherent drive and dissipation and showed that in the 
weak-pumping limit a simple perturbative-recursive ap- 
proach can be used to analytically compute the system 
density matrix elements as well as correlation functions. 
Such procedure can in principle be generalized to CQED 
arrays with more than two sites. The analytical solution 
was shown to be highly accurate over several decades of 
system parameters by comparing to the numerical solu- 
tion of the master equation for the system density matrix. 

Our results show that the zero-delay second-order cor- 
relation function provides an effective probe that can 
discriminate between coherent and strongly correlated 
regimes of a two-site CQED system. An interesting ex- 
tension of the present work is the generalization of the 
techniques and concepts used here to larger multi-site 
systems. 
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Appendix A: Density matrix elements 

We provide below the analytic solution of Eqs. (I7|-(l9| in the small F/^p limit, where we can truncate the polaritonic 

Hilbert space to Ntot = 2 and assume that the vacuum state to have approximately unit occupancy, i.e. pi^i ~ 1. 
The rest of the density matrix elements in steady state (see caption to Fig. [2] for the labelling of the basis states) is 
given below: 

Pi,i - 1 , (Al) 

Pl.2 = P2 1 = '^[JPl,^ + F{V2pi,i + Pl.l + Pl.5 - P2.2 - P3,2)] + 2%/2p2,4 + 2p3,5 , (A2) 

Pi, 3 = P3,l = ^[-^Pl.Z + ^(Pl,5 + \/2pi,6 + Pis - P2,3 " PZ,z)] + 2p2,5 + 2^2^3,6 , (A3) 

Pi,4 - /04,i = :^;^2iI7 [v^"^Pi,5 + F{V2pia - p2A - Pza)] , (A4) 

Pl.5 = P1.1 = ^ W^J{PlA + Plfi) + F{pi^3 + pi,2 - P2,5 - P3,5)] , (A5) 

Plfi = P*6,l = 7p-21J7 [^'^Pl,5 + F{V2pi^3 - P2.6 - PSfi)] , (A6) 

P2,2 = :^[JiP2,3 - Paa) + F{\/2p2A + P2.1 + P2,5 - Pl.2 - V^P4,2 " P5,2)] + 2p4,4 + P5.5 , (A7) 

/02,3 = P3,2 = ^[J{P2,2 - P3.3) + F{V2p2S + P2,l + P2,5 - Pl.3 - ^2^4,3 " P5,3)] + 2p4.4 + P5,5 , (A8) 

P3,3 = :^ [J{P3,2 - P2,3) + ■F'(\/2p3,6 + ^3,1 + P3,5 " Pl.3 " V^P6,3 " P5,3)] + 2p4.4 + P5,5 , (A9) 

P2A = P4,2 = 37pl21[/ ['^(^^^'5 " P3,4) + F{y/2p2,2 ~ Pl,4 " %/2p4,4 - P5,4)] , (AlO) 

P2.5 = P1^2 = 1^p^J{-s/^P2A + \/2p2,6 - P3,5) + ^(P2,3 + P2,2 ~ Pi, 5 - \/2p4,5 " Pb,b)\ , (All) 

P2,6 = PS,2 = a^^^air/ t-^^^^^^S - P3,6) + F{y/2p2,3 - Pi, 6 - %/2p4,6 - P5,6)] , (Al2) 

P3,4 = P4,3 = ^7p-2ir/ ['^(^^3,5 " P2.4) + -F(V2p3,2 " P5,4 ~ PlA ^ ^P6,4)] , (A13) 

P3,5 = P5,3 = 1^^[J{V2P3A + V2p3fi - P2,b) + ^(P3,3 + P3,2 " P5,5 " Pi, 5 " V2p6,5)] , (A14) 

P3,6 = P5,3 == f7p-21[/ ['^(^P3-5 ~ ^'2,6) + -P'(\/2p3,3 " P5,6 " %/2p6,6 ^ Pl,6)] , (A15) 

P4,4 = 2^ [\/2 J(P4.5 - P5,4) + \/2F(p4,2 - P2,4)] , (A16) 

P4,5 = P5,4 = 2-7P+21C/ [%/2 J(P4,4 + P4,6 - PS^s) + -F(P4,3 + P4,2 - ^/^P2,b)\ , (Al7) 

P4,6 = P6,4 = 2^ [\/2 J(P4,5 - P5,6) + \/2F(p4,3 ~ P2s)] , (A18) 

P5,5 = 2^ [V2 J(p5,4 + p5,6 - P6,5 " P4,5) + -F'(P5,3 + P5,2 - P3,5 " P2,5)] , (Al9) 

P5,6 = P6,5 = 27p-21J7 [^'^(P5,5 ^ P6,6 - P4,6) + ^(%/2p5,3 - P3,6 - P2,6)] , (A20) 

P6,6 = 2^ [V2 J(P6,5 - P5,6) + \/2F(p6,3 " P3,6)] , (A21) 

We simplify the equations above as follows. First we consider the equations with a first order dependence on F/"fp, 
and we neglect all the higher order terms. We get 

Pi,2 = p; 1 = — (Jpi.3 + ^Pi,i) , (A22) 

7p 



and 



This yields 



Pi,3 = p5 1 = — (Jpi.2 + Fpi,i) . (A23) 

1p 



Pi,2 = Pi,3 = P2,i = P3,i = -VTT2jT^-^Pi,i • (A24) 



.(1+21J) 

. ^ Ip ' 

+ (2J)2 - 

Next we insert these expressions into the remaining equations to obtain equations at a higher order in F/jp. Thus, 



^3 



the elements of p with a second order dependence on F/^p give the foUowing set of closed equations 

P2,2 = Y^[J{P2,?. - P3.2) + F{p2,i ~ Pl,2)] == ^F{p2s - Pi, 2) (A25) 

P2,3 == ::^[-'"(P2,2 - P3,3) + F{p2,l - Pl,z)] = -^F{p2,l - /Ol,3) (A26) 

P3,2 = ::^[-'"(P3,3 - P2,2) + -F(P3,1 " Pl^)] = ^^(/'S,! - Pl,2) (A27) 

P3,3 = ^^[J{ps,2 - P2.3) + F{p3^i - pi,3)] = j;F{p3.i - Pi.s) , (A28) 
from which wc can calculate a solution in terms of pi 1 as 

P2,2 = P2,3 = P3,2 = P3,3 = i + fily i^)^ Phi ■ (A29) 

Now we consider the equations of order (F/jp)^ 

PiA = Pli = ^^Pi.5 + ^^Pi.2 (A30) 

Pl,6 = Pll = :^P1,5 + ^^Pl.3 (A31) 

Pi,5 = Pli = 2^21^^1,4 + 2i^pi,2 . (A32) 

Again, we can solve the set of coupled equations reported above isolating the explicit dependence on /Oi.i, which gives 
the solutions 

-2y2(l+2M)2 ^ 

piA = Pifi = Pli = Pes = [T,^(^^_2iC/)+4j2]'[i+(|^)2]^ Phi . (A33) 

n - n* - 41J+2(7p-21C/) ^(^+^) f2^ ( aoa\ 

Ph5 - P5S - -4J2+7p(7^-21C/) 1+(^)2 ^ ^1^1 • l^^'^^'' 

We analyse the other set of equations of order (F/'jp)^ 

'°2,4 = l^^i^i^/ [>/(V2p2,5 - P2.6) + F{V2p2,2 - Ph^] , (A35) 

P^'^ - i^^i^it/ [^(^/2P2,5 - P3,6) + ^(V2p2,3 - Pl,6)] , (A36) 

P^'^ = l^^izir/ [^(^A3,5 - P2,4) + ^(V2p3,2 - Pl,4)] , (A37) 

P^^e = i^^i^ir/ [^(\/2p3,5 - P2,6) + ^(V2p3,3 - Pi,6)] , (A38) 

P2,5 - i :^ [ J(\/2/52,4 + V2p2fi - P3.5) + F(p2.,3 + P2a - Ph^)] , (A39) 

P3,5 = 14 [•^(V^P3,4 + V2p33 - P2.5) + i^(P3,3 + P3,2 - Pl,5)] , (A40) 

and we get a solution depending on the elements P2,2: Pi,4i Pi, 5; which in turn depend on pi 1 as shown above. The 
solutions for p2A ^nd p2,5 read 

„ _„ _„ _„ - n* - n* - n* - n* - i F (37p+2i./)(V2p2,2-pi,4)+2V21.7(2p2,2-pi,5) /A.iN 

P2,4 - P3,6 ~ P2,6 - P3,4 - P4,2 ^ P6,3 - P6,2 - P4,3 ^ ^^ (|7p-2i[/)(37p+2i J)+3i7p./+6J2 ' ^^^^) 

_ - n* - n* - i F 4V21J(^/2p2,2-Pl,4) + 2[(j7p-21i7) + lJ](2p2.2-pl,5) /AA9^ 

P2,5 - P3,5 - P5,2 - P5,3 - ^^ (|7p-2iC/)(37p+2i.7)+3i7p./+6J2 ' ^^^^^ 

We finally consider the terms depending on (F/jp)^, 

P4,5 = P5,4 = 27P+21C/ [V2J{2p4,4 - Ps^s) + ^(2p4,2 - \/2p2,5)] , (A43) 

P5,5 = ^[V2J{p5A - Ph5) + Fip5a - P2,5)] , (A44) 

P4,4 = 2^(72 J(p4,5 - P5,4) + \/2(p4,2 - P2,4) , (A45) 

P4,6 = 2^(\/2J(p4,5 - Pbfi) + \/2(p4,3 - P2,6) , (A46) 

P6,4 = 2^ (\/2 J(P6,5 - P5,4) + \/2(p6,2 - P3,4) , (A47) 

P6,6 = 2^ (\/2 J(P6,5 - P5,6) + ^2(^6,3 ^ P3,6) ■ (A48) 



from these closed set of equations we get a solution for ^4^4 as 



P4,4 -— P6,6 — P4,6 — P6A ~ 

+3(P2,5)8 J(2 J + [/) + 3ff(p2,5)87p J] , 



-3?(p2,4)8y27p>/ - 



(A49) 



where SR and 3 indicate real and imaginary parts of the respective elements. The explicit analytic expression for P4.4 
can be calculated from those of the elements P2.2, Pi, 4, Pi, 5, P2,4, P2,5, which have been obtained before. We notice 
that the recursive procedure described can be properly generalized to calculate the relevant density matrix elements 
also for generic multi-site CQED systems. 

Appendix B: Numerical solution 

To solve Eqs. (|7])-(|9| we use the finite-size Fock-state matrix representation of all the operators. For any given set 
of model parameters, the steady state density matrix can be obtained by numerically searching for the eigenvector 
|p))ss corresponding to the eigenvalue A^s = of the linear operator equation 



L\p)) = \\p)) 



(Bl) 



In the latter, jp)) is the density operator mapped into vectorial form, and L is the linear matrix corresponding to 
the Liouvillian superoperator on the right-hand side of the master equation. If it exists, as it is always the case for 
the parameters considered, the steady state solution is unique.^ After recasting the vector |p))ss in matrix form, 
the relevant observable quantities can be calculated as (O) ss — Tr\Opss\- In this work, we kept up to 6 photons 
per cavity in the basis, which is sufficient for convergence due to the weak driving conditions. Steady state results 
obtained in this way have been successfully compared to the ones obtained from a full time evolution of Eq. ([7]) as a 
further check. , 
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